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Abstract 

The evolutionary force of recombination is lacking in asexually repro- 
ducing populations. As a consequence, the population can suffer an ir- 
reversible accumulation of deleterious mutations, a phenomenon known 
as Muller's ratchet. We formulate discrete and continuous time versions 
of Muller's ratchet. Inspired by Ifaigh's (1978) analysis of a dynami- 
cal system which arises in the limit of large populations, we identify the 
parameter 7 — N\/{Ns ■ log(A''A)) as most important for the speed of 
accumulation of deleterious mutations. Here N is population size, s is the 
selection coefficient and A is the deleterious mutation rate. For large parts 
of the parameter range, measuring time in units of size A'^, deleterious mu- 
tations accumulate according to a power law in NX with exponent 7 if 
7 > 0.5. For 7 < 0.5 mutations cannot accumulate. We obtain diffusion 
approximations for three different parameter regimes, depending on the 
speed of the ratchet. Our approximations shed new light on analyses of 
Stephan et al. (1993) and Gordo & Charlesworth (2000). The heuristics 
leading to the approximations are supported by simulations. 

1 Introduction 

Muller's ratchet is a mechanism that has been suggested as an explanation for 
the evolution of sex (Maynard Smith 1978). The idea is simple; in an asexually 
reproducing population chromosomes are passed down as indivisible blocks and 
so the number of deleterious mutations accumulated along any ancestral line in 
the population can only increase. When everyone in the current 'best' class has 
accumulated at least one additional bad mutation then the minimum mutational 
load in the population increases: the ratchet clicks. In a sexually reproducing 
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population this is no longer the case; because of recombination between parental 
genomes a parent carrying a high mutational load can have offspring with fewer 
deleterious mutations. The high cost of sexual reproduction is thus offset by 
the benefits of inhibiting the ratchet. Equally the ratchet provides a possible 
explanation for the degeneration of Y chromosomes in sexual organisms (e.g. 
Charlcsworth 1978, 1996; Rice 1994; Charlcsworth & Charlcsworth 1997, 1998). 
However, in order to assess its real biological importance one should establish 
under what circumstances MuUer's ratchet will have an evolutionary effect. In 
particular, how many generations will it take for an asexually reproducing pop- 
ulation to lose its current best class? In other words, what is the rate of the 
ratchet? 

In spite of the substantial literature devoted to the ratchet (see Loewe 2006 
for an extensive bibliography), even in the simplest mathematical models a 
closed form expression for the rate remains elusive. Instead various approxima- 
tions have been proposed which fit well with simulations for particular parameter 
regimes. The analysis presented here unifies these approximations into a single 
framework and provides a more detailed mathematical understanding of their 
regions of validity. 

The simplest mathematical model of the ratchet was formulated in the pi- 
oneering work of Haigh (1978). Consider an asexual population of constant 
size A''. The population evolves according to classical Wright-Fisher dynamics. 
Thus each of the N individuals in the (t-l- l)st generation independently chooses 
a parent from the individuals in the tth generation. The probability that an 
individual which has accumulated k mutations is selected is proportional to its 
relative fitness, (1 — s)^ . The number of mutations carried by the offspring is 
then k -\- J where J is an (independent) Poisson random variable with mean A. 

Haigh identifies no = Ne~^^^ as an approximation (at large times) for the 
total number of individuals carrying the least number of mutations and finds 
numerical evidence of a linear relationship between uq and the average time 
between clicks of the ratchet, at least for 'intermediate' values of uq (which he 
quantifies as tiq > 1 and less than 25, say). On the other hand, for increasing 
values of no Stephan et al. (1993) note the increasing importance of s for the rate 
of the ratchet. The simulations of Gordo & Charlcsworth (2000) also suggest 
that for no fixed at a large value the ratchet can run at very different rates. 
They focus on parameter ranges that may be the most relevant to the problem 
of the degeneration of large non-recombining portions of chromosomes. 

In our approach we use diffusion approximations to identify another param- 
eter as being an important factor in determining the rate of the ratchet. We 
define 



' ' Ns \og{NX) ' ^ ' ' 

Notice that no = N{NX)~'^. In these parameters one can reinterpret Haigh's 
empirical results as saying that if we measure time in units of size N, then the 
rate of the ratchet follows a power law in A^A. In fact our main observation in 
this note is that for a substantial portion of parameter space (which wc shall 
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quantify a little more precisely later) we have the following 

Rule of Thumb. The rate of the ratchet is of the order N^-^X^ forj e (1/2, 1), 
whereas it is exponentially slow in (NX)^^'^ for 7 < 1/2. 

There are two novelties here. First, the abrupt change in behaviour at 
7=1/2 and second the power law interpretation of the rate for 7 G (1/2, 1). 
As an appetiser, Figure [T] illustrates that this behaviour really is reflected in 
simulated data; see also Jj6l 

The rule of thumb breaks down for two scenarios: first, if 7 > 1 then for 
large NX we have no < 1 and so our arguments, which are based on diffusion 
approximations for the size of the best class, will break down. This parameter 
regime, which leads to very frequent clicks of the ratchet, was studied by Gessler 
(1995). Second, if iVA is too large then we see the transition from exponentially 
rare clicks to frequent clicks takes place at larger values of 7. 

The rest of this note is laid out as follows. In ^ we review the work of 
Haigh (1978). Whereas Haigh's work focuses on discrete time dynamics, in fJ3]we 
write down instead a (continuous time) Fleming- Viot diffusion approximation 
to the model whose behaviour captures the dynamics of large populations when 
s and A are small. We then pass in 0to the infinite population (deterministic) 
limit. This system can be solved exactly and by, in the spirit of Haigh, using this 
deterministic system to estimate the behaviour of the 'bulk' of the population we 
obtain in fj5]a one dimensional diffusion which approximates the size of the best 
class in our Fleming- Viot system. The drift in this one-dimensional diffusion 
will take one of three forms depending upon whether the ratchet is clicking 
rarely, at a moderate rate or frequently per N generations (but always rarely 
per generation). Performing a scaling of this diffusion allows us to predict the 
relationship between the parameters NX and Ns of the biological model and 
the value of 7 at which we can expect to see the phase transition from rare 
clicking to power law behaviour of the rate of the ratchet. In f|6]we compare 
our predictions to simulated data, and in fj7]wc discuss the connection between 
our findings and previous work of Stephan ct al. (1993), Stephan & Kim (2002) 
and Gordo & Charlesworth (2000). 

2 The discrete ratchet — Haigh's approach 

The population dynamics described in the introduction can be reformulated 
mathematically as follows. Let iV be a fixed natural number (the population 
size), A > (the mutation parameter) and s G (0, 1) (the selection coefficient). 
The population is described by a stochastic process taking values in P(No), the 
simplex of probability weights on Nq. Suppose that x(<) = {xk{t))k=o,i,... S 
7'(No) is the vector of type frequencies (or frequency profile) in the tth gen- 
eration (so for example Nxk(t) individuals in the population carry exactly k 
mutations). Let H be an No-valued random variable with F[H = k] propor- 
tional to (1 — s)''xk{t), let J be a Poisson(A)-random variable independent of 
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(A) 



N=10' 




(B) 




Figure 1: (A) We plot the rate of the ratchet against 7, where time is measured 
in units of N generations. As predicted by our rule of thumb we see a sharp 
change of behaviour around 7 = 0.5. (B) We see the power-law behaviour of 
the rate for various values of 7. The solid lines are given by simulation of a 
Wright-Fisher model. The dashed lines fit the prediction that the time between 
clicks is a constant times N{NX)~'^ for 7 > 0.5. Note that this breaks down for 
7 > 1. 
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H and let Ki, be independent copies of iJ + J. Then the random type 
frequencies in the next generation arc 

Xk{t + l) = ^#{i:K\ = k}. (2.1) 

Wc shall refer to this as the ratchet dynamics in discrete time. 

First consider what happens as ^ oo. By the law of large numbers. (|2.ip 
results in the deterministic dynamics 

x(f + l) :=E,(t)[X(i + l)]. (2.2) 

An important property of the dynamics (|2.2p is that vectors of Poisson weights 
are mapped to vectors of Poisson weights. To sec this, note that when N oo 
the right hand side of (|2.2[1 is just the law of the random variable H+J. If x{t) = 
Poisson(a), the Poisson distribution with mean a, then H is Poisson («(! — s))- 
distributed and consequently x(i + 1) is the law of a Poisson(a(l — s) + A) 
random variable. We shall see in 31] that the same is true of the continuous time 
analogue of (|2.2p and indeed wc show there that for every initial condition with 
xq > the solution to the continuous time equation converges to the stationary 
point 

TT := Poisson((?) 

as t — > oo, where 




Haigh's analysis of the finite population model focusses on the number of 
individuals in the best class. Let us write k* for the number of mutations carried 
by individuals in this class. In a finite population, k* will increase with time, 
but the profile of frequencies relative to fc*, {^a;+a:»}, forms a recurrent Markov 
chain. We set 

Y :— {Yk)k=os... ■— {Xk'+k)k=o,i,.. 

and observe that since fitness is always measured relative to the mean fitness 
in the population, between clicks of the ratchet the equation for the dynamics 
of Y is precisely the same as that for X when Xq > 0. Suppose that after 
t generations there are no{t) ~ Nyo{t) individuals in the best class. Then the 
probability of sampling a parent from this class and not acquiring any additional 
mutations is 

Po{t) iyo{t)/Wit))e-\ 

where 

oo 

w{t) = J2y^m-sy■ (2.3) 

i=0 

Thus, given yoit), the size of the best class in the next generation has a binomial 
distribution with N trials and success probability po {t) , and so the evolution of 
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the best class is determined by W{t), the mean fitness of the population. We 
shall see this property of Yq reflected in the diffusion approximation of ^ 

A principal assumption of Haigh's analysis is that immediately before a click, 
the individuals of the current best class have all been distributed upon the other 
classes, in proportion to their Poisson weights. Thus immediately after a click 
he takes the type frequencies (relative to the new best class) to be 

n := -^(^1,^2, . . .) = f^^"'' T^"'' ■ (2.4) 

The time until the next click is then subdivided into two phases. During the first 
phase the deterministic dynamical system decays exponentially fast towards its 
Poisson equilibrium, swamping the randomness arising from the finite popula- 
tion size. At the time that he proposes for the end of the first phase the size of 
the best class is approximately I.Gttq. The mean fitness of the population has 
decreased by an amount which can also be readily estimated and, combining 
this with a Poisson approximation to the binomial distribution, Haigh proposes 
that (at least initially) during the second (longer) phase the size of the best class 
should be approximated by a Galton- Watson branching process with a Poisson 
offspring distribution. 

Haigh's original proposal was that since the mean fitness of the population 
(and consequently the mean number of offspring in the Galton Watson process) 
changes only slowly during the second phase it should be taken to be constant 
throughout that phase. Later refinements have modified Haigh's approach in 
two key ways. First they have worked with a diffusion approximation so that 
the Galton- Watson process is replaced by a Feller diffusion and second, instead 
of taking a constant drift, they look for a good approximation of the mean 
fitness given the size of the best class, resulting in a Feller diffusion with logistic 
growth. Our aim in the rest of this paper is to unify these approximations in a 
single mathematical framework and discuss them in the light of simulations. A 
crucial building block will be the following extension of (|2.4p , which we call the 
Poisson profile approximation (or PPA) of Y based on Yq: 

H(ro) := (Fo, ^^(^1,712,...)). (2.5) 

V 1 - TTo / 

As a first step, we now turn to a diffusion approximation for the full ratchet 
dynamics (|2.ip . 

3 The Fleming- Viot diffusion 

For large N and small A and s, the following stochastic dynamics on P(Mo) 
in continuous time captures the conditional expectation and variance of the 
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discrete dynamics (|2.ip : 
dXk= [Y. ~ k)XjXk + X{Xk-i - Xk) \dt + Y, JjJ^i^k dWjk, (3.1) 

fc = 0,l,2,..., 

where X-i :~ 0, and (Wjk)j>k is an array of independent standard Wiener 
processes with Wkj '■= —Wjk- This is, of course, just the infinite-dimensional 
version of the standard muhi-dimensional Wright-Fisher diffusion. Existence of 
a process solving p.ip can be established using a diffusion limit of the discrete 
dynamics of The coefficient s{j — k) is the fitness difference between type k 
and type j, A(Xj,„i — X^.) is the flow into and out of class k due to mutation 
and the diffusion coefficients -^XjXk reflect the covariances due to multinomial 
sampling. 

Remark 3.1. Often when one passes to a Fleming- Viot diffusion approxima- 
tion, one measures time in units of size N and correspondingly the parameters 
s and A appear as Ns and NX. Here we have not rescalcd time, hence the factor 
of jj in the noise and the unsealed parameters s and A in the equations. □ 

Writing 

3 

for the flrst moment of X, (|3.ip translates into 

dXk = (.s(Mi(X) - fc) - X)Xk + XXu-i)dt + Jlx^kdWjk (3.2) 

In exactly the same way as in our discrete stochastic system, writing k* for the 
number of mutations carried by individuals in the fittest class, one would like 
to think of the population as a travelling wave with profile Yk ~ Xi^^k' ■ Notice 
in particular that 



dYo = (sA/i(Y) - X)Yodt + y lro(l - Yo)dWo, 

where Wq is a standard Wiener process. Thus, just as in Haigh's setting, the 
frequency of the best class is determined by the mean fitness of the population. 
Substituting into (jS.ip one can obtain a stochastic equation for Mi, 



dMi = (A - sM2)dt + dG 
where M2 = J2jU " ^^^'^ the martingale G has quadratic variation 

d{G) = ^M2dt. 
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Thus the speed of the wave is determined by the variance of the profile. Simi- 
larly, 

= {-jjM2 + (A - sKh))dt + dH 

where M3 = Y^jij - MifXj and the martingale H has quadratic variation 
d{H) = ;^(^-'^4 — M2)dt with denoting the fourth centred moment and so 
on. 

These equations for the centred moments are entirely analogous to those 
obtained by Higgs & Woodcock (1995) except that in the Fleming- Viot setting 
they are exact. As pointed out there, Biirger (1991) obtained similar equations 
to study the evolution of polygenic traits. The difficulty in using these equations 
to study the rate of the ratchet is of course that they are not closed: the equation 
for Mk involves Mk+i and so on. Moreover, there is no obvious approximating 
closed system. By contrast, the infinite population limit, in which the noise is 
absent, turns out to have a closed solution. 



4 The infinite population limit 

The continuous time analogue of (j2.2p is the deterministic dynamical system 

dxk ^ ({s{Mi{yi) - k) ^ \)xk + \xk-i^dt, fc = 0,l,2,... (4.1) 

where X-i =0, obtained by letting ^ 00 in our Fleming- Viot diffusion (|3.2[) . 
Our goal in this section is to solve this system of equations. Note that Maia 
et al. (2003) have obtained a complete solution of the corresponding discrete 
system following (|2.2p . 

As we shall see in Proposition 14. 11 the stationary points of the system are 
exactly the same as for (|2.2p . that is x = tt and all its right shifts (tt^^^. )fe=o.i,... , 
fc* = 0, 1, 2, ... . Since the Poisson distribution can be characterised as the only 
distribution on Nq with all cumulants equal, it is natural to transform (|4.ip into 
a system of equations for the cumulants, k^,, fc = 1,2,..., of the vector x. The 
cumulants are defined by the relation 

logE-^^-^'-E«^^- (4-2) 

fe=0 k=l 

We assume > and set 

kq := -logXQ. (4.3) 

Proposition 4.1. For Kk,k = 0, 1, 2, . . . as in (|4.2p and (|4.3p the system (|4.ip 
is equivalent to 



-SKk+l + A, 



fc = 0,1,2,. 
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Setting k := (kq, ki, . . .) this system is solved by 

K = BK{Of + -{l-e-'')l, S = (6,y),.,=o,i,..., h,^\ 0"-)! 

s 10 otherwise. 

(4.4) 



In particular, 



and 

' A 



Ki(<) = ^fca.,(t)= - log^Xfe(0)e-«^- 

fe=0 ^ fe=0 



l-e""*). (4.6) 

s 

i=st 

Remark 4.2. If x(0) is a Poisson(yu) distribution then substituting into (|4.4p 
we see that x(i) is a Poisson distribution with parameter X/s + e^**(/i — A/s). 
In other words, just as for the discrete dynamical system considered by Haigh, 
vectors of Poisson weights are mapped to vectors of Poisson weights. In particu- 
lar TT := Poisson(A/s) is once again a stationary point of the system. Moreover, 
this proposition shows that for any vector x(0) with a;o(0) > 0, the solution 
converges to this stationary point. The corresponding convergence result in 
the discrete case is established in Maia et al. (2003). More generally, if k* is 
the smallest value of k for which Xk{0) > then the solution will converge to 

(7''/£-fc*)fc=0,l,2,...- 

Proof of Proposition \4-I\ Using (|4.ip we have 



A 



oo oo 



Xke 



d 



fc=0 

Thus, (1121) gives 

fc=0 fe=0 

Comparing coefficients in the last equation we obtain 

Kfe = -SKfc+l + A, fc:=0, 1, 
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This linear system can readily be solved. We write 
D := ((5i+ij)ij=o,i,2,..., 



so that 



Since 



= -sDk^ + Al^. (4.7) 



i — '- — J > I 



otherwise, 



the linear system (|4.7p is solved by 

kUV = e-^'^'KlO f + A / e-°''^ldu = e-^^^'KfO)^ + -(1 - e-"*)!^. 

Jo s 

□ 

Remark 4.3. With the initial condition x(0) :— tt given by p.4p . equations 
(|i3)) and (HH]) become 

-oW = e-^^--^ (4.8) 

and 

^i{t)-9~l+ (4.9) 

At time 

r:^'^, (4.10) 
s 

we have Xo{t) = i-e-i ~ LBttq. Comparing with SjH we see that in our 
continuous time setting r is precisely the counterpart of the time proposed by 
Haigh as the end of 'phase one'. 

In SjSlour prediction for A/i(Y) given Yq will require the value of Mi(y(T)) 
for y solving (|4.ip when started from a Poisson profile approximation. This is 
the purpose of the next proposition. 

Proposition 4.4. For i/q E (0, 1), let y{t) be the solution of (|4.ip with the 
initial state y(0) := n(yo) defined in (j2.5p . and let r be Haigh's relaxation time 
defined in ((jlO)l . Then for A>0 with 77 := 6'^-^ 



M.(y(A.)).. + ^(l-^) 
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Proof. Since 



^ TTfcg-^^ = exp ( - 9(1 - e-«)) = ttI 



we have 

oo 

k=0 

and 



fe=0 



^ ^ =yo + - 7ro(e -1 



= 2/0 + 7ro(e' - 1) 

^=sAt i — TTo 



d ^ _ 



k=0 



^=sAr 1 — TTq 



noe'r]. 



Using the solution (|4.5p and (|4.6p and j/o(0) = yo 



yoiAr) =2/0 



TToe'' 



2/0 



Mi(y(Ar)) 



2/o + S^'ro(e'7-l) 

7roe''(l - TTq) 

yo(l - TToe") + TTO (e" - 1)' 

- V 



yo + ^Me^'-l) 

TTo - 2/0 



'/ 



From ([4lT|) . 

and thus 
T^o - 2/0 



2/0 



yo(l-7roe'') + 7ro(e'?-l)' 
2/o(^T)7ro(e''-l) 



(4.11) 



(4.12) 



7roe"(l - TTo) - zjo(^t)(1 - ttoc") 
7roe''(7ro - 2/o(^t))(1 - ttq) 



7roe''(l - TTo) - 2/o(^t)(1 - TToe'?) ' 
yo{l - TToe") + 7ro(e" - 1) = 7ro(p/' - 1)- 



7roe''(l - TTo) 



7roe''(l - TTo) - 2/o(Ar)(l - TToe") ' 
Plugging the last two equations into (|4.12p we find 



Mi(y(Ar)) 



e '' - 1 



2/0 (^r) 



□ 



5 ONE DIMENSIONAL DIFFUSION APPROXIMATIONS 12 



p 
in 

O 
CO 

0.00 0.01 0.02 0.03 0.04 
Yo 

Figure 2: Using simulations (see also Q we plot {Yo,Mi). There is a good fit 
to a linear relationship between Yq and Mi. Note that 7 = 0.6 in the figure. 

5 One dimensional diffusion approximations 

Recall from that in our Fleming- Viot model the frequency Yq of the best 
class follows 



dYo^ (^sMi(Y)-X^Yodt+^^Yo{l~Yo)dWo, (5.1) 

where Wq is a standard Wiener process. The system of equations p.2p is too 
complex for us to be able to find an explicit expression for Mi ( Y) , which depends 
on the whole vector Y of class sizes. Instead we seek a good approximation of Mi 
given Yq. Substituting this into equation (|5.ip will then yield a one-dimensional 
diffusion which we use as an approximation for the size of the best class. Of 
course this assumption of a functional dependence between Yq and Mi is a 
weakness of the one-dimensional diffusion approximation, but simulations show 
that there is a substantial correlation between Yq and Mi, see, for example. 
Figure [H 

To understand our approach to finding a map Yq i— s- Afi, recall as a first 







N=10^ 
X=0.1 
s=0.024 






Yo=7ro 
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step Haigh's approximation that immediately after a click of the ratchet the 
profile has the form (|2.4p . The reasoning is as follows. Deviations of Y from a 
Poisson profile can only be due to the randomness arising from resampling in 
a finite population. Since resampling has no tendency to increase or decrease 
the frequency of a given class, the average profile immediately after a click of 
the ratchet is approximated by the state where ttq is distributed evenly over all 
other classes according to their equilibrium frequencies. During his short 'phase 
one', Haigh then allows this profile to 'relax' through the action of the discrete 
dynamical system (|2.2|) and it is the mean fitness in the population after this 
short relaxation time which determines the behaviour of the best class during 
'phase two'. 

A natural next step in extending this argument is to suppose that also in 
between click times the resampling distributes the mass ttq — Yq evenly on all 
other classes. In other words, given Yq, approximate the state of the system by 
Y = n(Fo) given by jM]). 

Of course in reality the dynamical system interacts with the resampling 
as it tries to restore the system to its Poisson equilibrium. If this restoring 
force is strong, just as in Haigh's approach one estimated mean fitness during 
phase two from the 'relaxed' profile, so here one should approximate the mean 
fitness Ml not from the PPA, but from states which arise by evolving the PPA 
using the dynamical system for a certain amount of time. We call the resulting 
states relaxed Poisson Profile approximations or RPPA. There are three different 
parameter regimes with which we shall be concerned. Each corresponds to a 
different value of rj in the functional relationship 

M,=e + -^(l-^). (5.2) 



e '' - 1 V TTo 

of Proposition 14.41 These can be distinguished as follows 



A small, Tj^e, Ml w (l->o), (5.3a) 

1 - TTo 

A^l, ?7 = 1, Ml w 6i + 0.58fl - — ), (5.3b) 



A large, w 0, Mik9+(i- —) (5.3c) 



The resulting maps Fo Mi are plotted in Figure [31 Observe that for consis- 
tency, Ml has to increase, on average, by 1 during one click of the ratchet. 

Finally, before we can apply our onc-dimcnsional diffusion approximation we 
must choose a starting value for Yq following equation (|5.ip . For A large, the 
system is already close to its new equilibrium at the time of a click and so we 
take Yq = ttq. 

For A = 1, at the time of the click we observe a state which has relaxed for 
time T from a state of the form tt from (|2.4[) . We computed in Remark 14.31 that 
such a state comes with Yq — I.Gttq. 



5 ONE DIMENSIONAL DIFFUSION APPROXIMATIONS 



14 




Figure 3: Since simulations show a strong correlation between the first moment 
Ml and Yq, we use (|5.3a|) - (|5.3cP to predict Mi from Yo depending on the model 
parameters. 



For small values of A, observe that the profile of the population immediately 
after a click is approximately tt from (|2.4p . Since tt is not a state of the form 
n(yo) the arguments that led to Proposition l4.4l do not apply. Instead we follow 
Haigh in dividing the time between clicks into two phases. Consider first 'phase 
one'. Recall from the dynamical system that 

dYo = [sMi - X)Yo dt. 

We write 

(sMi - X)Yo = c (tto - Yo), 

where c (like Yq and Mi) depends on r = 6e^^* . Starting in Y(0) = tt we have 
from and gH) 

6le-"* « r 

>o(t) = 



1 - exp(-6'e-«*) 1 - e- 



A/i(t) =9-1 + ' =9-1 

cxp{9e ^'^j — 1 



We compute 



c (Ml -9)Yo _ 1 - _ r(l 
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It can be checked that this expression hes between 1 and 1.25 for all r > 
which suggests that the size of the best class in the initial phase after a click is 
reasonably described by the dynamics 



dYo = s (tto - yo) dt + ^ ^Yo dWo (5.4) 

started from jj^^^- We allow Yq to evolve according to equation (|5.4p until 
it reaches I.Gtto, say, and then use our estimate of Mi from equation ()5.3p to 
estimate the evolution of Yq during the (longer) 'phase two'. 

We assume that states of the ratchet are RPPAs, i.e., Poisson profile ap- 
proximations (|2.5p which are relaxed for time At, where t = ^ log 6, which 
leads to the functional relationship (|5.2p . Consequently, we suggest that (|5.ip 
is approximated by the 'mean reversion' dynamics 



dYo = (1 - ^)Yodt + sj ]^Yo dWo, (5.5) 

with rj = 9^^^ , where we have used a Feller noise instead of the Wright-Fisher 
term in ()5.ip . In other words, Yo is a Feller branching diffusion with logistic 
growth. 

Using the three regimes from (j5.3p . we have the approximations 



A small, dYo = A(7ro - ^0)^0^^ + \/ ^YodW, (5.6a) 



A = l, dYo = 0.b%s[l-^\Yodt + J —YodWo, (5.6b) 

* TTo / \ N 



A large, dYo = s{l- ^^Yodt + ^Yo dWo, (5.6c) 

(where in the first equation we have used that j-rj- ~ 1 + ttq and that Yotto is 
negligible). 

An equation similar to (|5.6b[) was found (by different means) by Stephan et 
al (1993) and further discussed in Gordo & Charlesworth (2000). Stephan and 
Kim (2002) analyse whether a prefactor of 0.5 or 0.6 in (|5.6bp fits better with 
simulated data. We discuss the relationship with these papers in detail in ^JT) 

The expected time to extinction of a diffusion following (|5.5[) is readily ob- 
tained from a Green function calculation similar to that in Lambert (2005). 
We refrain from doing this here, but instead use a scaling argument to identify 
parameter ranges for which the ratchet clicks and to give evidence for the rule 
of thumb formulated in the introduction. 



Consider the rescaling 



Z{t) = —Yo {Nnot) . 
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For A small equation (|5.6ap becomes 

dZ ^ NXir^il - Z)Zdt + y/ZdW 

= {N\f~^'<{l - Z)Zdt + \fZdW. (5.7) 

For ^ = 1 on the other hand we obtain from (|5.6bp 

dZ = 0.587Vs7ro(l - Z)Zdt + ^/^dW 

= {).b&—-^——{N\f-'^{l-Z)Zdt + yfZdW. (5.8) 
7log(A' A) 

For A large we obtain fr'om (|5.6cp the same equation without the factor of 0.58. 

From this rescaling we see that the equation that applies for small A, i.e. 
(|5.7p . is strongly mean reverting for 7 < 1/2. Recall that the choice of small A 
is appropriate when the ratchet is clicking frequently and so this indicates that 
frequent clicking simply will not happen for 7 < 1/2. To indicate the boimdary 
between rare and moderate clicking, equation (|5.8p is much more relevant than 
equation (|5.7p . At first sight, equation (|5.8p looks strongly mean reverting for all 
7 < 1, which would seem to suggest that the ratchet will click only exponentially 
slowly in (A^A)^~^. However, the closer 7 is to one, the larger the value of A^A 
we must take for this asymptotic regime to provide a good approximation. For 
example, in the table below we describe parameter combinations for which the 
coefhcient in front of the mean reversion term in equation (|5.8p is at least five. 
We see that for 7 < 1/2 this coefficient is large for most of the reasonable values 
of iVA, whereas for 7 > 1/2 it is rather small over a large range of A^A. 



7 


0.3 


0.4 


0.5 


0.55 


0.6 


0.7 


0.8 


0.9 


A^A > 


20 


102 


9 • 102 


4-103 


2 • 10"^ 


4-106 


2 - 10" 


8 - 1026 



Thus, for example, if 7 = 0.7 we require A' A to be of the order of IO6 in order 
for the strong mean reversion of equation (|5.8p to be evident. This is not a 
value of A^A which will be observed in practice. Indeed, as a 'rule of thumb', for 
biologically realistic parameter values, we should expect the transition from no 
clicks to a moderate rate of clicks to take place at around 7 ~ 0.5. 



6 Simulations 

We have argued that the one-dimensional diffusions (|5.6p approximate the fre- 
quency in the best class and from this deduced the rule of thumb from 21 In 
this section we use simulations to test the validity of our arguments. 

For a population following the dynamics (|2.ip . the {t + l)st generation is 
formed by multinomial sampling of A^ individuals with weights 

^^(^) = ^ w{t) ' 7!' ^'-'^ 
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where W{t) is the average fitness in the tth generation from (|2.3p and it is this 
Wright-Fisher model which was implemented in the simulations. 

To supplement the numerical results of Figurc[T]wc provide simulation results 
for the average time between clicks (where time is measured in units of N 
generations) for fixed N and 7 and varying A; see Figure S) Note that, for fixed 
7 in equation (jl.ip . s is increasing with A . We carry out simulations using 
a population size of TV = 10^ and A varying from 10~^ to 1. For 7 — 0.5 we 
observe that the power law behaviour breaks down already for TV A lO'^ and 
the diffusion (j5.6bp predicts the clicking of the ratchet sufficiently well. For 
increasing 7, the power law breaks down only for larger values of A'^A. For 7 = 
0.7, in our simulations we only observe the power law behaviour but conjecture 
that for larger values of NX the power law would break down; compare with the 
table above. 

For a finer analysis of which of the equations (|5.6|) works best, we study the 
resulting Green functions numerically; see Figure [S] In particular, we record 
the relative time spent in some cIYq in simulations and compare this quantity 
to the numerically integrated, normalised Green functions given through the 
diffusions (|5.6ap and (|5.6bp . (We do not consider (|5.6cp because it only gives 
an approximation if the ratchet clicks rarely.) We see in (A) that for 7 = 0.5 
not only does (j5.6bp produce better estimates for the average time between 
clicks (Figured]) but also for the time spent around some point yo. However, 
for 7 = 0.9, clicks are more frequent and we expect (|5.6ap to provide a better 
approximation. Indeed, although both (|5.6aP and (|5.6bp predict the power law 
behaviour, as (B) shows, the first equation produces better estimates for the 
relative amount of time spent in some (IYq. 

To support our claim that the states observed are relaxed Poisson Profile 
Approximations we use a phase-plane analysis; see Figure [6l At any point in 
time of a simulation, values for Yq and Mi can be observed. The resulting 
plots indicate that we can distinguish the three parameter regimes introduced 
in ^ In the case of rapid clicking of the ratchet (so that the states we observe 
have not relaxed a lot and thus are approximately of the form Il{yo)) we see 
in (A), (B) that the system is driven by the restoring force to Mi < 9. The 
reason is that Mi is small at click times and these are frequent. However, the 
slope of the line relating Yq and Mi is low, as predicted by (|5.3ap . (We used 
(1 — Yq) /(I — ttq) ~ 1 — Yq + ttq in the plot here.) For the case A = 1 the system 
spends some time near Yq ^ and thus the ratchet clicks, but not frequently. 
So, the dynamical system restores states partly to equilibrium and we see that 
the slope given in (|5.3bp gives the most reasonable prediction in (C), (D), (E). 
For rare clicking, i.e. A large, the dynamical system has even more time and 
(F) shows that the slope is as predicted by (|5.3cp . 

Our prediction that we observe profiles which are well approximated by a 
relaxed PPA applies especially well at click times. We check this numerically by 
observing the frequency of the (new) best class at click times; see Figure [Tj For 
small 7 the ratchet clicks rarely and the system has some time to relax to its 
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Figure 4: The power law behaviour of the rate of the ratchet with respect to 7 is 
vahd for a large portion of the parameter space. (A) For 7 = 0.5 clicks become 
rare and the power law does not apply for A^A > 10^. (B), (C) For 7 = 0.55 and 
7 = 0.6, we have to explore a larger portion of the parameter space in order to 
see that the power law does not apply any more. (D), (E), (F) For 7 > 0.7 we 
never observe a deviation from the power law. For every plot, wc used = 10^ 
and simulations ran for 5 • 10^ (7 = 0.5: 2 • 10'') generations for each value of 
A^A. 
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Figure 5: We compare the plots for the occupation density of Yq from the 
simulations with theoretical curves corresponding to the Green functions for 
the cases of small A and A = 1 in (|5.6p . (A) If clicks are rare, A = 1 produces 
better results than small A. (B) If clicks are frequent, the simulated densities of 
Yq are better approximated by small A. Every plot is based on the simulation 
of 5 • 10^ generations. 

new equilibrium even before the click of the ratchet. As a consequence, we see 
that the frequency of the (new) best class at the time of the click is already close 
to ttq. However, if 7 is large and the ratchet clicks frequently, the dynamical 
system has no time before the click to relax the system to the new equilibrium. 
Therefore, we observe that the frequency of the new best class is close to tti. 

7 Discussion 

Haigh (1978) was the first to attempt a rigorous mathematical analysis of the 
ideas of MuUer (1964). However, in spite of the apparent simplicity of Haigh's 
mathematical formulation of the model, the exact rate of Muller's ratchet re- 
mains elusive. In this note, we have developed arguments in the spirit of 
Haigh (1978), Stephan et al. (1993) and Gordo & Charlesworth (2000) to give 
approximations for this rate. 

Haigh gave the empirical formula 

4^7ro + 71og6l+|-20 

for the average time between clicks of the ratchet (where time is measured 
in generations). A quantitative understanding of the rate was first obtained 
by Stephan et al. (1993) using diffusion approximations and later extended by 
Gordo & Charlesworth (2000). Both obtain the diffusion (j5.6bp as the main 
equation giving a valid approximation for the frequency path of the best class. 
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Figure 6: There are three regimes for the relationship between Yq and Mi, as 
given in (j5.3p . If cHcks are frequent, at least the slope of the relationship between 
Yq and Mi fits roughly to equation (|5.3ap . If clicks occur reasonably often, 
(j5.3b[) gives a good approximation. If clicks are rare, (|5.3c|) gives a reasonable 
prediction. The plots show simulations for different values of 7. The dashed 
horizontal and vertical lines are Mi = 9 and Yq = ttq, respectively. For every 
plot we used N = 10"', A = 0.1 and simulations ran for 10^ generations. 
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Figure 7: Our heuristic that observed states come from a RPPA apphcs in 
particular at chck times. (A) For smaU 7, i.e., rare chcking, the frequency of 
the best class is already close to ttq while (B) it is close to tti for larger 7, i.e., 
frequent clicking. Every plot is based on the simulation of 5 • 10^ generations. 



The reasoning leading to (j5.6bp in these papers is twofold. Stephan et 
al. (1993) and Stephan & Kim (2002) argue that although fitness decreases 
by se~^ during one 'cycle' of the discrete ratchet model from ^ (in which the 
system advances from one Poisson equilibrium to the next), at the actual click 
time only a fraction of the fitness has been lost. They suggest kse~^ for k = 0.5 
01 k = 0.6 as the loss of fitness at click times. In other words they predict the 
functional relationship Mi{Yo) discussed in ^by linear interpolation between 
Mi(7ro) = 6 and Mi(0) = 9 + ke~^ Ki + k; compare with Figure [3l On the 
other hand, Gordo & Charlesworth (2000) use a calculation of Haigh which tells 
us that if the dynamical system (j2.ip is started in tt from p.4p . then at the 
end of phase one (corresponding in the continuous setting, as we observed in 
Remark 1131 to time i log 6*) we have sMi « 1 - e"-^(l + 0.42s). This leads to 
the approximation Mi(1.67ro) = — 0.42 and again interpolating linearly using 
Mi(7ro) = 9 gives (|5.6bp . 

Simulations show that (|5.6bp provides a good approximation to the rate of 
the ratchet for a wide range of parameters; see e.g. Stephan and Kim (2002). 
The novelty in our work is that we derive (|5.6b[) explicitly from the dynamical 
system. In particular, we do not use a linear approximation, but instead derive 
a functional linear relationship in Proposition 14.41 In addition, we clarify the 
role of the two different phases suggested by Haigh. As simulations show, since 
phase one is fast, it is already complete at the time when phase two starts. 
Therefore, in practice, we observe states that are relaxed PPAs. 

The drawback of our analysis is that we cannot give good arguments for the 
choice oi A = 1 in (j5.3bp and (|5.6bp . However, note that the choice of A = 1 is 
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essential to obtain the prefactor of 0.58 in (j5.3bp . E.g., if = 10, the choice of 
A = 0.5 leads to a prefactor of 0.13 while A = 2 leads to the prefactor of 0.95, 
neither of which fits with simulated data; see Figure [HI 

We obtain two more diffusion approximations, which are valid in the cases of 
frequent and rare clicking, respectively. In practice, both play little role in the 
prediction of the rate of the ratchet. For fast clicking, (|5.6b|) shows the same 
power law behaviour as (j5.6ap and rare clicks are never observed in simulations. 

Of course from a biological perspective our mathematical model is very naive. 
In particular, it is unnatural to suppose that each new mutation confers the 
same selective disadvantage and, indeed, not all mutations will be deleterious. 
Moreover, if one is to argue that MuUer's ratchet explains the evolution of sex, 
then one has to quantify the effect of recombination. Such questions provide a 
rich, but challenging, mathematical playground. 

Acknov^rledgements We have discussed MuUer's ratchet with many different 
people. Wc arc especially indebted to Ellen Baakc, Nick Barton, Matthias 
Birkncr, Charles Cuthbcrtson, Don Dawson, Wolfgang Stcphan, Jay Taylor and 
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